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The local structure of a solid-on-solid (SOS) interface in a two-dimensional kinetic Ising ferromag- 
net with single-spin-flip Glauber dynamics, which is driven far from equilibrium by an applied field, 
is studied by an analytic mean-field, nonlinear-response theory [P. A. Rikvold and M. Kolesik, J. 
Stat. Phys. 100, 377 (2000)] and by dynamic Monte Carlo simulations. The probability density of 
the height of an individual step in the surface is obtained, both analytically and by simulation. The 
width of the probability density is found to increase dramatically with the magnitude of the applied 
field, with close agreement between the theoretical predictions and the simulation results. Excellent 
. agreement between theory and simulations is also found for the field-dependence and anisotropy of 

the interface velocity. The joint distribution of nearest-neighbor step heights is obtained by simula- 
tion. It shows increasing correlations with increasing field, similar to the skewness observed in other 
examples of growing surfaces. 
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The motion of surfaces and interfaces plays a central role in many scientific and technological disciplines. In 
particular, the dynamics of interfaces such as phase- and grain boundaries in solid materials [Jy and domain walls 
in magnets Q and ferroelectrics heavily influence both dynamic and static materials properties. Among interfaces 
characteristic of two-dimensional systems are steps on crystal surfaces || , domain walls in thin magnetic and dielectric 
films H , and boundaries between different types of vegetation such as savanna and rainforest Q . 
■ An enormous amount of work in recent years has been devoted to the dynamics and structure of moving and growing 
interfaces However, despite the fact that many important interface properties, such as mobility and catalytic 

and chemical activity, are largely determined by the microscopic interface structure, most of this effort has been 
concentrated on large-scale scaling properties. Although the detailed atomistic mechanisms by which interfaces move 
are often not known, useful understanding can be obtained from stochastic models in which the motion occurs through 
random nucleation and migration of local topological features such as kinks and steps . It is therefore important 
to gain better insight into how, for different stochastic dynamics, the driving force (such as an applied magnetic or 
electric field, a chemical-potential gradient, or the amount of rainfall in the case of models of vegetation distribution) 
may alter the microscopic structure of the interface, thereby leading to a highly nonlinear velocity response. 

In a previous paper [Q we introduced a dynamic mean-field approximation for the microstructure of an interface 
in a two-dimensional kinetic Ising ferromagnet with a single-spin-flip Glauber dynamic j^, S|, driven by an applied 
field |n| |Tl[ . This model is directly applicable to many magnetic and ferroelectric systems and other cases where the 
interface dynamics are not inhibited by coupling to a conserved field [ 12| , Based on the resulting local interface 
structure, we obtained a nonlinear-response approximation for the steady-state propagation velocity, which was shown 
to be in good agreement with dynamic Monte Carlo (MC) simulations for a wide range of fields and temperatures. 
However, since the approximation was based on the Burton-Cabrera-Frank (BCF) unrestricted solid-on-solid (SOS) 
model (lj| , the overhangs and bubbles in the Ising interface were handled in an uncontrolled way. Here we therefore 
consider the performance of our approximation for the unrestricted SOS model, so that overhangs and bubbles are 
absent at all times by definition of the model. In particular, we obtain the surface velocity under the Glauber dynamic 
as a function of applied field and temperature, as well as its anisotropy for tilt angles between 0° and 45°. 
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In a recent paper [[L5| we showed that the microscopic interface structure, and thus the mobility, can depend 
dramatically on the details of the dynamics. The most significant difference is between dynamics in which the 
transition probabilities of the local spins factorize into one part that depends only on the change in interaction energy 
due to the transition and one that depends only on the change in the field energy (soft dynamics [^6| ), and dynamics 
that cannot be factorized in this way (hard dynamics Mj). For soft dynamics, the interface structure for all values 
of the field remains the same as in equilibrium at zero field |jl5|| , so that the mobility can be calculated exactly at the 
level of linear response. In the present paper we concentrate on the standard Glauber dynamic |8|, ^] , which belongs 
to the class of hard dynamics and thus leads to more complicated and interesting behavior. 

Both the driven Ising and SOS surfaces belong to the Kardar-Parisi-Zhang (KPZ) dynamic universality class || |l7| , 
in which the macroscopic, stationary distribution for flat, moving interfaces is Gaussian, corresponding to a random 
walk with independent increments. Nevertheless, the step heights in several discrete models in this class are correlated 
at short distances [l9| . In the mean-field approximation developed here, these short-range correlations are ignored. 
The resulting discrepancies, which are minor, are elucidated by comparison with MC simulations. 

The remainder of this paper is organized as follows. In Sec. II we introduce the SOS interface model and derive 
a linear-response approximation for its velocity in a nonzero external field. In Sec. Ill we develop a mean-field 
approximation for the time evolution of the single-step probability density function (pdf), as well as for its stationary 
form. The latter enables us to extend the approximation for the interface structure and velocity to a nonlinear- 
response level. The analytical approximations are compared with dynamic MC simulations in Sec. [V. In Sec. IV A 



we numerically solve the mean-field equation of motion for the single-step pdf and compare the resulting values 
of the time-dependent average interface step height with MC si mulati ons. In Sec. IV B we compare the simulated 
stationary single-step pdfs with the theoretical predictions. In Sec. IV C we compare the simulated stationary interface 
velo city w ith the theoretical predictions for various values of applied field, temperature, and interface tilt angle. In 
Sec. IV D we compare simulations and analytical predictions for the detailed stationary interface structure, including 
the asymmetry of the simulated interface at nonzero fields. A summary and conclusions are found in Sec. 



II. MODEL AND DYNAMICS 



The original BCF SOS model considers an interface in a lattice gas or S = 1/2 Ising system on a square lattice 
of unit lattice constant as a series of integer- valued steps 5(x), parallel to the y axis. The interface height is thus 
a single-valued integer function of the x-coordinate, with steps at integer values of x. The configuration energy is 
given by the nearest-neighbor Ising Hamiltonian with anisotropic, ferromagnetic interactions J x and J y in the x and 
y direction, respectively: 

7~L = ~ ^ S Xy y ( J x S x +l,y + JyS x ,y+l + H) . (1) 

Here the two states at lattice site (x, y) are denoted by s XtV = ±1, and runs over all sites. The quantity H is the 
applied field, and the interface is introduced by fixing s XiV — +1 and —1 for large negative and positive y, respectively. 
Without loss of generality we take H > 0, such that the interface on average moves in the positive y direction. 

A single-spin-flip (nonconservative) dynamic which satisfies detailed balance, such as the Metropolis or Glauber 
algorithms j^, ^j, ensures the approach to equilibrium, which in this case is a uniformly positive phase with the 
interface pushed off to positive infinity. Such algorithms are defined by a single-spin transition probability, W(s XtV — > 
—s x ,y) = W(/3AE, (3U). Here (3 is the inverse of the temperature T (Boltzmann's constant is taken as unity), AE is 
the energy change corresponding to a successful spin flip, and the optional parameter U is an energy barrier between 
the two states that enters into Arrhenius-type stochastic dynamics [^0| and other dynamics that include a transition 
state The detailed-balance condition is expressed as W(f3AE, (3U) /W(—f3AE, (3U) = e~' 3AB , where the right- 
hand side is independent of U. (In the case of soft dynamics, the detailed-balance condition is satisfied independently 
for the two parts of W.) In order to preserve the SOS configuration at all times, flips are allowed only at sites which 
have exactly one broken bond in the y direction. 

With the Ising Hamiltonian there are only a finite number of different values of AE, and the spins can therefore be 
divided into classes |ll], , labeled by the spin value s and the number of broken bonds between the spin and its 

nearest neighbors in the x and y direction, j and fc, respectively. The ten spin classes consistent with the SOS model 
are denoted jks with j £ {0, 1, 2} and fee {0,1}. They are listed in Table |. 

In this paper, as in Ref. M, we use the standard discrete-time Glauber dynamic with the transition probability 
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Time is measured in MC steps per spin (MCSS). 

In the BCF SOS model the heights of the individual steps are assumed to be statistically independent and identically 
distributed. The step-height probability density function (pdf) is given by the interaction energy corresponding to 
the \S(x)\ broken J^-bonds between spins in the columns centered at (x — 1/2) and (x + 1/2) as 



p[6(x)] = Z^r 1 ^^! e 7W)5(a;) 



(3) 



The Boltzmann factor X — e~ 2 P Jx determines the width of the pdf, and 7(0) is a Lagrange multiplier which maintains 
the mean step height at an x- independent value, (S(x)) = tan0, where is the average angle between the interface 
and the x axis. The partition function is 



Z(0) = ]T X^e 



1-X 1 



with 7(0) given by 



5— — OC 



1 + X 2 ) tan. 



1 - 2X cosh 7 (0) +X 2 



(4) 



(1-X 2 



tan" 1 



AX' 



1/2 



2X{1 + tan0) 

(see details in Ref. 0). For = 0, 7(0) = 0, and the partition function simplifies to 

Z(0) = (l + X)/(1-X). 



(5) 



(0) 



The mean spin-class populations, (n(jks)), are all obtained from the product of the independent pdfs for S(x) and 
S(x+l). Symmetry of p[<5(a;)] under the transformation (x,4>, S) — > (— x, — 0, — S) ensures that (n(jk— )) = (n(jk+)} 
for all 7 and k. Numerical results illustrating the breakdown of this up/down symmetry for large \H\ are discussed 
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As discussed in Ref. J7J , calculation of the individual class populations is straightforward but somewhat 
tedious, especially for nonzero 0. The final results are summarized in Table [0|. 

Whenever a spin flips from —1 to +1, the corresponding column of the interface advances by one lattice constant 
in the y direction. Conversely, the column recedes by one lattice constant when a spin flips from +1 to —1. The 
corresponding energy changes are given in the third column in Table |. Since the spin-class populations on both sides 
of the interface are equal in this approximation, the contribution to the mean velocity in the y direction from sites in 
the classes jk— and jk+ becomes 



(v y (jk)} = W (/3AE(jk-),j3U) - W (0AE(jk+),pU) 

The results corresponding to the hard Glauber transition probabilities used here, Eq. 
of Table ||. The mean propagation velocity perpendicular to the interface becomes 



(7) 



are listed in the last column 



(t>^(n(iks))(v y {jk)) 



(8) 



where the sum runs over the classes included in the tables. While the general result is cumbersome if written out in 
detail, the special case of = leads to a compact formula: 



(v ± (T,H,0)) 



tanh(/3if) 
(1 + AO 2 



2X 



1 + X 2 



sinh(2f3J m ) 
cosh(f3H) 



It was shown in Ref. that Eq. (||) reduces to the results for the single-step O, 
growth Jl^, |2^, [27|] models at low temperatures for strong and weak fields, respectively. 



(9) 
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III. NONLINEAR- RESPONSE 

With X = e~ 2l3J " , the results in Table || correspond to a linear-response approximation. In Ref. Q] we developed 
a mean- field approximation leading to a field-dependent X(T,H) 1 based on a detailed-balance argument for the 
stationary state. Here we show that this detailed-balance relation follows naturally from a dynamic mean-field 
approximation for the equation of motion for the single-step pdf during the approach to the stationary state. 
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We denote the total transition probability for the height of the single step at x to change from 8(x) to 8(x) ± 1 as 
W[<5(x) — > S(x) ±1]. In terms of these transition probabilities, the equation of motion for the single-step pdf, p[S(x), t] 
becomes 

dp[ y = -p[8(x),t}{W[S(x)^8(x)-l] 

+W[5(x) -» «(jb) + 1]} 
+p[*(ar) + 1]W[<%) + 1 -► <5(x)] 

+p[5(s) - l]W[S(x) - 1 -» *(as)] , (10) 

where the coupling to the joint multistep probability density is hidden in the single-step transition rates W. 

To obtain an approximation for W[8 — > S ± 1], we employ the same mean-field assumption of independent steps as 
in equilibrium. For S(x) > 1 to increase to 8(x) + 1, cither the spin in front of the interface at £+1/2 can flip from 
— 1 to +1, or the spin behind the interface at a;— 1/2 can flip from +1 to —1. In each of these cases, AE can have two 
different values, depending on the value of 5(x+l) and 8(x— 1), respectively. The same argument holds also for the 
reverse transition, 8(x) + 1 — * 6(x). The energy changes and corresponding conditions on 8(x+l) and 8(x— 1), which 
are shown in Fig. [j], yield 

W[5^S±l,t] = ~{\W(-2H) + W(+2H)]IL±(t) 

+ {W(-2H + 4J X ) + W(+2H + 4J X )][1 - II±(t)]} 

VV[*±1 -►*,<] = i{[W(-2if) + W(+2if)][l-n T (i)] 

+ [W(-2iT-4J x )+W(+2if-4J x )]II T (t)} , (11) 

where the upper signs refer to 8 > +1 and the lower signs to 8 < —1, 11+ (t) = 5Z£f +1 p[<M], andll_(t) = 5Zs^lif>[<M]- 
For simplicity, we here write the single-site transition rates W(f3AE, j3U) as W(AE). 
In the stationary limit, Eqs. (10) and (g) lead to the detailed-balance condition 



o±lW _p[8±l] _W[5->5±1] 



X(T, me^W = i^-r^ = —-77 - £ , (12) 

v ' ' p[8] W[8±l ->■£]' v ' 



where the upper and lower signs have the same interpretation as in Eq. (jllj). Using Eqs. (|3|) and (|4|), we get the 
stationary values for 11+ and n_ , 



n. 



Xe^ (1- le^(«) 



1 -X 
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n = (13) 



1-X 2 



which, when inserted together with Eq. (O) in Eq. (12|), yield a self-consistency equation for X. The self-consistency 
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equation reduces to a linear equation for X 2 , and wit 
rates W, the solution takes the form 



1 the help of the detailed-balance condition for the spin transition 



1/2 



X(T H) = f e~ 2 " HW (- 2H ~ 4J -) + ^ H W(+2H - 4J X ) ■ 

V ' ; \ W(-2if - 4J X ) + W(+2H - 4 J x ) J { ' 

which is independent of 7(0). (This solution is the same as the one obtained in Ref. [Q, but the derivation given in 
that paper did not explicitly show that all dependence on "f(4>) cancels out.) 

Equation ( [Til) shows that X(T,H) depends on the specific dynamic, except for H — 0, where it reduces to its 
equilibrium value, X(T, 0) = e _2/3Ja! , for all dynamics that satisfy detailed balance. A situation in which the in- 
dependence in Eq. (jlj) cancels out, is that of the soft dynamics discussed in Sec. I. (The barrier energy U can 
be contained in one or the other of the factors.) In Ref. [|l5| we demonstrated that a soft dynamic yields an SOS 
interface that is identical to the equilibrium SOS interface at if = and the same temperature, regardless of the 
value of H . Neither the Glauber dynamic used here, nor the equally common Metropolis dynamic with transition 
probability Wm{sx,v — * — s x, y ) = mm [l, e~ /3A£; ] ||, g], satisfies this factorization condition. Such hard dynamics lead 
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to a nontrivial field dependence in X. Inserting the Glauber dynamic defined by Eq. (||) into Eq. (|Tj), we explicitly 
get 

2BJ f e 2 ^ cosh(2/3H) + e" 2 ^! 1/2 , , 

All the results for the spin-class populations of the zero-field equilibrium interface, which are listed in Table 0, 
can now be applied to obtain a nonlinear-response approximation for the steady-state propagation velocity of flat, 
driven interfaces with hard dynamics. This simply requires replacing the zero-field X = e~ 2 ^ ,7j; used in the linear- 
response approximation by the field-dependent X(T,H), obtained from Eq. ( |l4| ) using the transition probabilities 
corresponding to the particular dynamic used. For soft dynamics the linear-response result with X = e~ 2l3J:c is exact 
©■ 

A physical reason for the marked difference between hard and soft dynamics is best seen by comparing concrete 
examples of dynamics in the two classes, such as the hard Glauber dynamics used here, Eq. (||), and the soft Glauber 
dynamics used in Ref. Jl5|| , 

W SG (s x , y - -s x>y ) = [1 + eP AE »] 1 ■ [1 + e^ E -'\ 1 , (16) 

in the case of a very strong field. In the hard case, the effect of the field completely dominates the transition rates, such 
that the rate is near unity for transitions that bring a spin parallel to the applied field, and near zero for transitions 
in the opposite direction, irrespective of the change in interaction energy. In the soft case, the probability of bringing 
a spin antiparallel to the field is also near zero, but the probabilities of different transitions bringing a spin parallel 
to the field differ according to the corresponding change in the interaction energy, as given by the second factor in 
Eq. @- 

In the next Section we show that this nonlinear-response approximation gives very good agreement with MC 
simulations of driven, flat SOS interfaces evolving under the hard Glauber dynamic for a wide range of fields and 
temperatures. 



IV. COMPARISON WITH MONTE CARLO SIMULATIONS 



We have compared the analytical estimates of step-height distributions, propagation velocities, and spin-class 
populations developed above with MC simulations of the same model for J x = J y = J. The details of our particular 
implementation of the discrete-time n-fold way rejection- free MC algorithm |2^] are the same as described in Ref. [Q, 
except that only transitions from the classes with one broken y-bond (k = 1) are allowed. By keeping only the 
interface sites in memory, the algorithm is not subject to any size restriction in the y direction, and simulations can 
be carried out for arbitrarily long times. 

The numerical results presented here are based on MC simulations at the two temperatures, T = 0.2T C and 0.6T C 
(T c = — 2J/ln(V2 — 1) is the critical temperature for the isotropic, square-lattice Ising model pg[]), with L x = 10 000 
and fixed <p between and 45°. In order to ensure stationarity we ran the simulation for 5 000 n-fold way updates per 
updatable spin (UPS) before taking any measurements (100 000 UPS for some of the strongest fields and largest values 
of 4> at T — 0. 2T e ). Exploratory simulations with both larger and smaller L x (up to 100 000) and "warm-up" times 
(see Sec. iVA) showed that the values used in the production runs were sufficient to ensure a stationary interface. 
Stationary class populations and interface velocities were averaged over 50000 UPS. In the stationary limit 1 UPS 
corresponds to between 2 MCSS for strong fields at both temperatures, and about 75 MCSS for H = at T — 0.2T C . 
Adequate statistics for one- and two-step pdfs was ensured by the large L x , ten times the value used in Ref. 0. 
Each data point took approximately nine hours on a DEC Alpha Unix workstation or four hours on a Pentium II PC 
running Linux. 



A. Approach to stationarity 

In order both to check the applicability of the mean-field approximation at early times, and to decide the approxi- 
mate time needed to reach the stationary state, we first studied the transient behavior of the average step height (|<5|) 
for cj) = at T = 0.2 J and 0.6J. 

In the dynamic MC simulation (\5\) was measured during a "time window" which was opened after a specified 
number of UPS, corresponding to a given approximate average evolution time t in MCSS. To obtain an optimum 
balance between time resolution and accuracy, the width of the window was varied from approximately one MCSS 
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at early times, to about ten MCSS at late times, and five independent runs were performed for each value of t. 
Standard errors for t and (|<5|) were estimated in the usual way from the spread in their measured values over the five 
realizations. 

For comparison with the MC simulations, we solved the mean-field equation of motion for the single-step pdf, 
Eq. ([ll]), numerically by a first-order iterative scheme with a time step of 10~ 4 MCSS (shorter time steps made no 
discernible difference). Both the simulations and the solution of the equation of motion were started from a sharp 
interface at t = 0. Results for fields between and 10 J are shown in Fig. [|. 

In general, we find overall qualitative agreement between the simulations and the equation of motion. For T = 0.2T C 
and \H\ < 2 J, both methods have reached a common stationary value by t = 10 000 MCSS, while for T = 0.6T C 
and \H\ < 3 J, stationarity is reached by t — 1000 MCSS. Our choice of 5 000 UPS as "warm-up time" in our 
studies of the stationary properties, corresponding to at least 10 000 MCSS, is thus well justified. However, there are 
significant quantitative differences between the simulation results and the solution of Eq. ( |l0| ) for early times. We 
believe this indicates that the mean-field assumption of statistically independent step heights is not well justified until 
the interface structure has been "randomized" through a sufficient number of updates. For the extremely strong field, 
H = 10 J, W(S — > S ± 1) pa W(S ± 1 — * 8) w 1/2, so that the evolution of the interface width is essentially diffusional 
and (\S\) oc t 1 / 2 . Both methods agree with this result, although the amplitude for the MC simulation is larger than 
predicted by the mean-field equation of motion. 



B. Single-step probability densities 



Single-step pdfs were obtained by MC simulation at T = 0.2T C and 0.6T C for 4> = and several values of H between 
and 3.0J. The simulation results for p[S] are shown in Fig. ||, together with the theoretical result, Eq. (|J) with 
X(T, H) from Eq. (|l5|). For both temperatures, the agreement is excellent in the whole range of S and H shown. 

A simple comparison between the analytical and simulation results is given in Fig. ^(a), which shows (\S\) vs H for 
4> = at T = 0.2T C and 0.6T C . The solid curves represent the theoretical result obtained by summation of Eq. (|^), 
(\S\) = 2Xj (l — AT 2 ), with X from Eq. (|l5|). There is excellent agreement between the theoretical field dependence 
and the MC data. Additional confirmation of the form of the single-step pdf, Eq. (||) , is obtained from the simulation 
results by calculating both directly by summation over the numerically obtained pdf and from the probability 
of zero step height as (|<5|) = {vW 1 ~ p[0]} /2. 

A slightly different way to check the agreement between the analytical predictions and the simulation results for 
the single-step pdf, is to compare X(T, H) as given by Eq. ( p"5| ) with the same quantity obtained from the simulations 
under the assumption that Eq. (|J) holds. From this equation for the pdf, using Z(0) from Eq. (0), it follows that 
X is given in terms of p[0] as X = {1 -p[0]}/{l + »[0]}, and in terms of (\S\) as X = y/1 + (\S\)- 2 - Both 
these MC estimates for X(T,H) are shown in Fig. [|(b). Again, the agreement is excellent. The slight deviations of 
the estimate based on (\S\) for large H are probably due to the fact that data were only recorded in separate bins for 
\5\ < 64, so that the calculated average becomes inaccurate whenever higher steps cannot be ignored. The estimate 
based on p[0] does not suffer from this problem. However, the slight discrepancy between both MC estimates on the 
one hand and Eq. (15) on the other, which is seen between H/J = 0.5 and 2 for T = 0.2T C , is probably a real effect. 



C. Interface velocities 

In this section we compare the simulated interface velocities with the analytical approximation, Eq. (||) . Figure |5| 
shows the normal velocity vs H for </> = 0. Included are both the linear-response approximation (i.e., X — e~ 2 P Jl ) 
and the nonlinear-response result with X(T,H) from Eq. (|l5|). Overall, there is excellent agreement between the 
MC results and the nonlinear-response theory, while the linear-response approximation seriously underestimates the 
velocity, especially at the lower temperature. As for the Ising model with hard Glauber dynamics studied in Ref. [Q, 
these results show that the latter approximation is clearly inadequate, and we include no further linear-response 
results in this paper. A slight disagreement between the simulations and the analytical predi ctions is seen at 0.2T C in 



the same range of field values as for X(T, H) obtained from the single-step pdf in Sec. |VB. 

The dependence of the normal velocity on the tilt angle <fi is shown in Fig. |^ for several values of H/J between 
0.1 and 3.0. At T = 0.2T C the anisotropy undergoes a gradual change from increasing with in agreement with the 
polynuclear growth model at small angles and the single-step model for larger angles at weak fields, to Eden- type 
inverse anisotropy [^9[ ||(], ||l| at strong fields [Fig. ^(a)]. At T = 0.6T C , on the other hand, inverse anisotropy is found 
for the stronger fields, growing gradually more pronounced with increasing H [Fig. ||(b)], while for the weakest fields 
studied the velocity is nonmonotonic in tan</) [Fig. ^|(c)]. The agreement between the simulations and the analytical 
results is excellent everywhere. 
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The temperature dependence of the normal interface velocity is shown in Fig. |7] for several values of H/J between 
0.2 and 3.0. The agreement between the simulations and the analytical results is excellent everywhere. 

D. Spin-class populations and skewness 

A closer look at the performance of the mean-field approximation for the interface structure is provided by the 
mean spin-class populations. The analytical predictions for the class populations are based on the assumption that 
different steps are statistically independent. A comparison of the simulation results with the analytical predictions 
therefore gives a way of testing this assumption. 

The six mean class populations, (n(01s)), (n(lls)), and (n(21s)) with s = ±1 are shown vs H in Fig. |8] for = 
and T — 0.2T C and 0.6T C . At both temperatures the analytical approximations follow the average of the populations 
for s = +1 and s = — 1 quite well, but at intermediate fields in particular, the populations in front of the surface 
(s = —1) and behind it (s = +1) are distinctly different. 

The skewness between the spin populations on the leading and trailing edges of the interface are a consequence of 
short-range correlations between neighboring steps, and it is quite commonly observed in driven interfaces. This is 
the case, even when the long-range correlations vanish as they do for interfaces in the KPZ dynamic universality class, 
to which the present model belongs. Such skewness was also observed in our study of the Ising model in Ref. (j, but 
in that case it was difficult to separate it from the effects of bubbles and overhangs. Skewness has also been observed 
in other SOS-type models, such as the body-centered SOS model studied by Neergaard and den Nijs Q and a model 
for the local time horizon in parallel MC simulations studied by Korniss et al. (However, no skewness is observed 
for the soft Glauber dynamic, a result which may be general for soft dynamics.) The correlations associated with the 
skewness generally lead to a broadening of protrusions on the leading edge ("hilltops"), while those on the trailing edge 
("valley bottoms") are sharpened pl| , or the other way around Jl9| . In terms of spin-class populations, the former 
corresponds to (n(21— )) > (n(21+)) and (n(ll+)) > (n(ll — )). The relative skewness can therefore be quantified by 
the two functions, 

_ <n(21-)) - <n(21+)) 
9 <n(21-)) + (n(21+)) ' 1 J 

introduced in Ref. Q], and 

(n(ll+)) - (n(ll-)) 
- („(ll+)) + <n(ll-)) ■ 1 ] 

These two skewness parameters are shown together in Fig. |^. The relative skewness is seen to be considerably stronger 
at the lower temperature. This temperature dependence is especially pronounced for p. 

Yet another way to visualize the skewness is to consider the joint two-step pdf, p[5(x),5(x + 1)]. Logarithmic 
contour plots of this quantity for different values of H are shown in Fig. [To] for = at T — 0.6T C . It is clearly 
seen how the contours change with H. For H—0 a symmetric diamond shape with equidistant contours indicates 
statistical independence with single-step pdfs given by Eq. (||). For stronger fields we find shapes that are elongated 
in the second quadrant [S(x) < 0, S(x + 1) > 0] and foreshortened in the fourth quadrant [S(x) > 0, 8(x + 1) < 0] . This 
shape indicates that large negative 8(x) tend to be followed by large positive 8(x + 1) (sharp valleys), while positive 
5(x) tend to be followed by smaller negative 8{x + 1) (rounded hilltops). We have not been successful in attempts to 
construct an analytical approximation which describes this evolution of the joint two-step pdf with increasing field. 

V. CONCLUSION 

In this paper we have considered in detail the microstructure of an unrestricted solid-on-solid (SOS) interface 
with Glauber dynamics, which is driven far from equilibrium by an applied field. The microstructure is of interest 
because it determines a number of interface properties, such as mobility and chemical reactivity. We adapted to 
this model a mean-field, nonlinear-response approximation previously developed for driven Ising interfaces without 
the SOS restriction |Q. In comparison to the Ising driven interface, which leaves bubbles of the unstable phase in 
its wake and exhibits "overhangs," the SOS interface is a relatively simple object. The absence of overhangs and of 
fluctuations in the stable and unstable phases (bubbles behind and in front of the interface) makes the SOS interface 
more suitable for description in terms of a mean-field type model. Moreover, unlike the Ising model, in which there 
are several effects that simultaneously contribute to the inaccuracy of the approximate treatment, the simpler SOS 
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structure makes it possible to identify the short-range correlations as the only significant factor causing deviations 
between the true interface behavior and the mean-field theory. 

To study the microstructure of the interface in detail, we investigated the interface velocity as a function of driving 
field, temperature, and angle relative to the lattice axes. We also studied the local shape of the interface in terms of 
the spin-class populations and the probability density for individual steps in the interface. In essentially all cases we 
found excellent agreement between our theoretical description of the stationary moving interface and the results of 
our dynamic MC simulations. 

The microstructure of the moving interface depends crucially on the details of the stochastic dynamics, and for 
the Glauber dynamic used here, the average height of a step in the interface was found to increase strongly with 
the applied field. Our theory predicts that this should be the case (with quantitative variations depending on the 
particular dynamic) for any dynamic in which the single-spin transition rates cannot be factorized into a part that 
depends only on interaction energies and another that depends only on the applied field (hard dynamics In 
contrast, for factorizing (soft) dynamics the interface structure should remain independent of the field. This was 
recently confirmed for the unrestricted SOS interface by MC simulations ]l5| . It is therefore important to use great 
caution in drawing conclusions about the microstructure of driven interfaces, based on dynamic MC simulations. For 
such conclusions to be valid, the dynamic must be chosen appropriately to the physical system of interest. The hard 
type of dynamics would appear to be particularly suited for certain interfaces in magnetic or dielectric systems, where 
the local order parameter is not conserved. 

One aspect of the interface dynamics not completely captured by our model is, of course, the short-range correlations. 
Namely, within the mean-field approximation used here, individual steps of the interface are assumed to be statistically 
independent. However, for increasing fields the interface undergoes a gradual breakdown of up/down symmetry. This 
is clearly seen in our simulations here, as well as in other examples of driven interfaces p8| jfjl. It would seem likely 
that one could construct a mean-field approximation at the two-spin level, which might be able to predict this skewness 
for hard dynamics, as well as its absence for soft dynamics. However, such a theory has not yet been developed. 
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TABLE I: The spin classes in the anisotropic square-lattice SOS model. The first column contains the class labels, jks. The 
second column contains the total field and interaction energy for a spin in each class, E(jks), relative to the energy of the state 
with all spins parallel and H — 0, Eo = —2(J X + J y ). The third column contains the change in the total system energy resulting 
from reversal of a spin from s to — s, AE(jks). In both E(jks) — Eo and AE(jks), the upper sign corresponds to s = —1, and 
the lower sign to s — +1. The first three classes (marked *) have nonzero populations in the SOS model, and flipping a spin 
in any of them preserves the SOS configuration. The other two classes also have nonzero populations in the SOS model, but 
flipping a spin in any of them would produce an overhang or a bubble and is therefore forbidden. 



Class, jks 


E(jks) - E 


AE(jks) 


01s * 


±H + 2Jy 


T2H + 4J X 


lis * 


±H + 2{J X + Jy) 


=F2# 


21s * 


±H + 2{2J X + J y ) 


T2H - 4J X 


10s 


±H + 2J X 


T2H + 4Jy 


20s 


±H + 4J X 


=F2# - 4(J X - Jy) 



TABLE II: The mean populations for the spin classes of the SOS interface, with the corresponding contributions to the 
interface velocity under the hard Glauber dynamic. The first column contains the class labels, jks. The second column 
contains the mean spin-class populations for general tilt angle 0, with cosh "/(</)) from Eq. (|^). The third column contains 
the spin-class populations for = 0. Using X = e~ 2l3J:c in these expressions yields the linear-response result in which the 
spin-class populations are evaluated for H = 0. Using X = X{T,H) from Eq. (jjj) with the transition probabilities of the 
particular dynamic used [here: Glauber, Eq. dTHI)], one gets the nonlinear-response approximation. The fourth column contains 
the contributions to the mean interface velocity in the y direction from spins in classes jk— and jk+, Eq. ((?]), using the 
SOS-preserving hard Glauber dynamic. 



Class, jks 


(n(jks)), general cf> 


{n(jks)),<t> = 


{Vy(jk)) 


01s * 


1-2X cosh ~/(4>)+X 2 


i 


tanh(/3H) 




(l+X) 2 


1 + 


- S inh(2.J /, >- 
o.;-ll< iH) 


2 


lis * 


2X\(1 + X 2 ) cosh -y(<t>) -2X] 


2X 


tanh (J3H) 


(1-X^Y 2 


(l+X) 2 


21s * 


X'- ! [1-2X cosIitOHX^] 


X 2 


tanh(/3H) 


(l-X^) 2 


(l+X) 2 


1 + 


r S inh(2.Jj, i" 


2 


10s 


2X 2 I 2cosh 2 y(<p)-l-2X cosh j(<f,) + X 2 


2X 2 (1+2X) 





1-X' 2 1 1-2X cosh j(<f,) + X 2 

X 2 [1-2X cosh7(,£) + X 2 ] 1 
(1-X 2 ) 2 ] 


(1-X 2 )(1+X) 2 


20s 


X' l [l-2X co5Yi~f(4>)+X 2 } 


X 4 





(1-X 2 ) 3 


(1-X 2 )(1+X) 2 
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6(-l)-l; 8(0)+l 5(0) >0 5(0)+l; S(+l)-l 



FIG. 1: Figure for calculating the single-step transition rates in Eq. ([ll|). Interface configurations are shown by bold line 
segments. Spins above (in front of) the interface equal —1, and spins below (behind) the interface equal +1. At the center is 
shown a step 5(0) > 1 (here shown as 5(0) = +1). A transition to 5(0) + 1 can be effected by flipping either of the two spins 
in the dashed boxes. The resulting configurations, which depend on the heights of the neighboring steps, are shown to the 
right and left in the figure. The corresponding energy changes and conditions on the neighboring step heights are given next 
to the arrows. The arrows pointing outward from the center of the figure correspond to the transition described above, while 
the arrows pointing toward the center correspond to the reverse transition, 5(0) + 1 — > 5(0). After Ref. [q|. 
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FIG. 2: The mean step height (\5\), shown vs time t (in MCSS) as predicted by the mean-field equation of motion for the 
single-step pdf, Eq. (jlO) (solid curves), and by dynamic MC simulations (crossed error bars indicating statistical standard 
errors in t and (|<5|)). From bottom to top, the results shown are for H/J — 0, 1, 2, 3, and 10 [MC results for H/J = 10 in part 
(b) only], (a) T = 0.2T C . (b) T = 0.6T C . 
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FIG. 3: MC (data points) and analytical (solid lines) results for the stationary single-step pdf, shown on logarithmic scale vs 
S. The fields are H/J — (filled squares), 0.5 (filled triangles), 1.0 (filled stars), 1.5 (filled diamonds), 2.0 (empty squares), 2.5 
(empty triangles), and 3.0 (empty stars), (a) T = 0.2T C . (b) T = 0.6T C . 
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FIG. 4: (a) Average stationary step height (\S\) vs H for (j>—0 at T=0.2T C and 0.6T C . The curves represent the theoretical 
result. The MC data were obtained directly by summation over the simulated single-step pdfs (filled symbols) and from the 
probability of zero step height (empty symbols). See text for details. Curve with filled circles and empty squares: 0.2T C . 
Curve with filled triangles and empty circles: 0.6T C . (b) The stationary pdf width parameter X(T,H) vs H: the analytical 
result Eq. (|l^) (curves) and estimates based on MC simulation results with 0=0 for p[0] (empty symbols) and for (\S\) (filled 
symbols). See text for details. Curves and symbols have the same interpretations as in (a). In this and all the following figures, 
the statistical uncertainty is much smaller than the symbol size. 



13 




0.5 1 1.5 2 2.5 3 



H/J 
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FIG. 5: The average stationary normal interface velocity {v±} vs H for <j> = 0. The MC results are shown as data points, the 
linear-response results as dashed curves, and the nonlinear-response results as solid curves, (a) T = 0.2T C . (b) T = 0.6T C . 
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FIG. 6: The average stationary normal interface velocity (v±) vs tan0 for several values of H. MC data are represented by 
data points, and analytical results by solid curves, (a) T = 0.2T C . (b) T = 0.6T C . In parts (a) and (b), the values of H/J are 
(from below to above): 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 [in (b) only], (c) H/J = 0.1 for T = 0.6T C , shown on a magnified 
scale to reveal the nonmonotonic dependence on tan (j> for this very weak field. 
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0.2 0.4 0.6 0.8 J 
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FIG. 7: The average stationary normal interface velocity (v±) vs T for <j> = 0. MC data are represented by data points, and 
analytical results by solid curves. From below to above, the values of H/J are 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0. 
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FIG. 8: Mean stationary class populations (n(jks)} vs H for (j> — 0. (a) T = 0.2T C . (b) T = 0.6T C . From top to bottom at 
the left edge of both parts, the classes are 01s, lis, and 21s with squares representing MC data for s = +1 and triangles for 
s = — 1. The analytic approximations are indicated by the solid curves. Note the different vertical scales in the two parts. 
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FIG. 9: The two relative skewness parameters p (triangles) and e (squares), defined in Eqs. ([I?]) and (|is{), respectively. The 
parameters are shown vs H for (j> — 0. (a) T = 0.2T C . (b) T — 0.6T C . Note the different vertical scales in the two parts. 
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FIG. 10: Contour plots of log 10 p[5(x),S(x + 1)] for = 0. (a) H/J = 0. (b) H/J = 1.0. (c) ff/J = 2.0. Note the different 
scales in the three parts. See discussion in the text. 
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